#include <stdio.h>
#include <stdlib.h>
#include <omp.h>
#include <math.h>


void tweakedScanBitOps(double *X, int n) {
	int i, j, d, logn, step, limit, mask, bit;
	//TYPE *c = (TYPE *)X;

	logn = (int)ceil(log2(n));
	bit = 1;
	mask = 0;

	//double t1 = omp_get_wtime();
	for (d = 1; d <= logn; ++d) {
		mask |= bit;  // increase mask
		step = bit << 1;
//		printf("mask = %X \n", mask);
//		printf("bit = %X \n", bit);

#pragma omp parallel for shared(X, mask, bit, n) private(j) schedule(guided) 
		for (i = mask; i < n; i += step) // has to have at least mask to enter if.
		{
		// select with lower bits
//			if ((i & mask) == mask) {
				// get the paired element.  j is always smaller than i because of the bit.
//				j = i ^ bit;
				// now do the math
				X[i] += X[i^bit];
//			} else {
//				printf("error \n" );
//			}
		}
		bit = step;
	}


	mask = (1 << logn) - 1;
	bit = 1 << (logn-1);
	
	for (d = logn-1; d > 0; --d) {
		mask ^= bit;  // keep lower bits only
		step = bit;
		bit >>= 1;   // use to find the half step merge target.
//		printf("mask = %X \n", mask);
//		printf("bit = %X \n", bit);
		limit = n - bit;
		
		// update i, i is the paired target.  j has the computed sum that
		// needs to be added to i. e.g. i = 3, and j = 2.  this updates value
		// of i.
#pragma omp parallel for shared(X, limit, n, step, bit, mask) private(j) schedule(guided) 
		for (i = mask; i < limit; i += step)  // has to have at least mask, and can't go over n-bit.
		{
			// select lower bits using mask
//			if ((i & mask) == mask) {
				// get the paired target
//				j = i + bit;
//				if (j < n) {
					X[i+bit] += X[i];
//				} else {
//					printf("ERROR: j = %d and limit = %d\n", j, limit);
//					printf("       i = %d and bit = %d (%X)\n", i, bit, bit);
//				}
//			} else {
//				printf("error \n" );
//			}
			
		}
	}

	//double t2 = omp_get_wtime();

	//printf("\n\t inside scan : time = %g", (t2-t1));
}
